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I. INTRODUCTION 


Setniclassical treatments of quantum systems are always of 
broad interests. Examples are nuclei, metallic clusters, cold 
atomic gases, neutron star crusts, etc. This is particularly 
useful for large systems which challenge the capacities of su¬ 
percomputers. In this context, the Thomas-Fermi approxima¬ 
tion (or Local Density Approximation) has been extended to 
higher orders to include as much as possible quantum cor¬ 
rections. When pairing correlations are included, the ex¬ 
tended Thomas-Fermi method for superconducting/superfluid 
systems becomes much more complicated than for normal 
fluid systems. In the past few years, only a very limited num¬ 
ber of works concerning the formulism of the second-order 
superfluid Thomas-Fermi approach have been given in the lit- 
erature|l2-[^ with either no or very limited practical calcula¬ 
tions Q. So further studies and applications are very desir¬ 
able. Furthermore, the density-dependent effective mass and 
the spin-orbit potential have not been considered yet, which 
are essential ingredients in, e.g., generalized nuclear density 
functionals such as the widely used Skyrme nuclear density 
functionals 01. Also in cold Fermi gases, the spin-orbit 
coupling |1 and the density-dependent effective mass 1 are 
currently very interesting. 

The Hartree-Fock-Bogoliubov (HFB) equations H, or 
Bogoliubov-de Gennes equations d, have been a framework 
for superfluid Fermi systems. The first S-expansion of the 
HFB density matrix is based on the Wigner Kirkwood trans¬ 
formation of the Bloch propagator ll, however, without prac¬ 
tical applications. Later, the /i-expansion of the Green’s func¬ 
tion of HFB solutions has been derived for superconduct¬ 
ing systems ll. More recently, a third paper based on the 
Green’s function method appeared with some applications to 
cold atomic systems 1. We also noticed that the coarse grain¬ 
ing treatment in superfluid Fermi gases M has a connection 
to the second-order Thomas-Fermi method. 

Our objective of the present work is to extend the for¬ 
mulism to include the density-dependent effective mass and 
the spin-orbit effect together with the second-order Thomas 
Fermi approximation for superfluid systems. Although such 


extensions d and even tensor interactions in have been 
achieved long time ago for non-superfluid systems, they are 
not yet elaborated for superfluid systems. In this work, we first 
derive the full correction terms based on the ^-expansion of 
the Green’s function corresponding to the HFB equations. In 
addition, we implement these second-order correction terms 
in numerical examinations and compare the results with fully 
self-consistent Skyrme HFB calculations of nuclei where the 
density-dependent effective mass and the spin-orbit potential 
are included. In the past the Thomas-Fermi approximation 
with Gogn y fo rce has already been applied to superfluid nu¬ 
clei in din. 

Recently, the adaptive multi-resolution 3D coordinate- 
space HFB method has been developed for complex superfluid 
systems but the computation of continuum states is still very 
costly S. The coordinate-space HFB calculations are very 
useful for describing weakly-bound systems and complex¬ 
shaped systems in, in contrast to the conventional HFB ap¬ 
proaches based on harmonic-oscillator basis expansions. A 
promising way to address large systems is to adopt the hy¬ 
brid HFB calculations m, i.e., using the Thomas-Fermi ap¬ 
proximation for high-energy quasi-continuum states and the 
discretized HFB solutions for low-energy states. The hybrid 
HFB has first been adopted with the zeroth-order Thomas- 
Fermi method for cold atomic systems dill. In this re¬ 
spect, large superfluid systems such as the trapped cold atoms 
with 10^“® particles 17^ |2^ and exotic neutron star crusts 
in large 3D cells id |22|] are numerically very challenging. 
In such inhomogeneous systems, the violation of zeroth-order 
Thomas-Fermi approximation related to finite-size effects and 
complex spatial topologies could be non-negligible Idl . 


II. FORMULATION 


The HFB e quat ion in the coordinate-space representation 
takes the form 1241] : 
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where h = /ihf — A; /ihf is the Hartree-Fock Hamiltonian; 
A is the Fermi energy (or chemical potential); A is the pairing 
potential; Uk and 14 are the upper and lower components of 
quasi-particle wave functions, respectively; is the quasi¬ 
particle energy. 

The theoretical derivation of the second-order superfluid 
Thomas-Fermi approximation starts with the solution of the 
HFB equation using the Green’s function method (or Gorkov 
equation) H,!!], 


■ / h{r) A(r) 

VA*(r) -h[r) 


G{r,r',(jj) = f^6{r — r'). (2) 


The Wigner transformation of the Green’s function form of 
the HFB equation can be written as; 


The various expansion terms of G can be obtained order by 
order. Compared to Ref.||2l, we include the spin-orbit poten¬ 
tial entailing more terms to appear in Eq.®. First, Go cor¬ 
responds to the zeroth-order solutions and does not change 
due to the spin-orbit potential. The matrix elements of Go are 
written as. 
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where the Gg^^^ and Gg^'*'^ correspond to the normal density 
and pairing density, respectively. The abbreviation D for the 
denominator in Eq.® is defined as D = — A^). 

Gi changes due to the spin-orbit potential and is written as 


A[(u;I - H(R, p)), G(R', p', cu)] = hi, (3) 

where H refers to the HEB Hamiltonian, and the operator A 
can be expanded in powers of h: 


Gusli* — -f^hso[{uj + hY - A^] 

1 ( 10 ) 

= ^h,o{2hA) 
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Similar to non-superfluid systems, the spin-orbit potential 
can be included as a higher-order term in terms of h in the 
Hamiltonian Q, 


Note that in this work we are treating systems in the case of 
real pairing potentials. The formulism with complex pairing 
potentials can be derived similarly but would be more compli¬ 
cated ||2] . 

G 2 is quite complicated 13, Ht] and the additional terms due 
to the spin-orbit potential are: 


H(R,p) = H,(R,p)+H,o(R,p), (5) 

in which and H^jo denote the central term and the spin-orbit 
term of the HEB Hamiltonian, respectively. The spin-orbit 
Hamiltonian is defined as 


Gi'i’= + '■)’-(" +S'.) A^l 


( 11 ) 
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The extended Thomas-Eermi method consists in expanding 
the Green’s function in Eq.® in powers of h , i.e., 

00 

G(R,p,a;) = ;i^G„(R,p,w)/i’^. ( 7 ) 

n—0 

The expansion of Eq.® in terms of H to the second-order 
can be written as, 

(cul — Hc)Go = I, (8a) 


Note that terms involving derivatives of the spin-orbit poten¬ 
tial have been omitted. Actually the terms in G 2 with hso do 
not contribute and only terms with do. 

Next, the contributions to normal density and pairing den¬ 
sity can be obtained by an appropriate integration of G 2 in the 
complex oj plane. Eor convenience of writing the expressions, 
we introduce a differential operator A as in Ref. H, 

^ = (Vr'^p-V p^R), (12) 

which is the operator in the Poisson bracket. Ea. (fT^ can be 
applied repeatedly and is related to Eq.® without expansion 
coefficients. The 2"‘^-order normal density contribution terms 
are given as. 
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and the 2''‘^-order pairing density, 


;52(R,p) = 0i(/iA2/i) + 02(/iA2A) + 03(AA2A) 
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(14) 

In Eg. (foi l and Eg.jTO. the coefficients rji and 9i have been 
given in Refs. (ll-IS, although they have different conven¬ 
tions. We refer to the Appendix for explicit expressions of 
these terms. Note that some terms vanish with the application 
of operators Aj, since we suppose a momentum-independent 
pairing potential. The additional contributions due to the spin- 
orbit potential are given in Eas. (fT3T l and (fT4l) explicitly, in 
which E denotes the HFB quasiparticle energies and takes the 
form as 


E = V(/ihf(R, p)-A)2+A(R, p)2. (15) 


where {Eq, E^) or (eg, £c) defines an energy interval in which 
the Thomas-Fermi approximation is used. 

The 2"‘^-order kinetic density correction in real spaces is 
obtained as. 


Mr) 



/•^° P2(R,p) P^E 
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(19) 


The 2'''^-order spin-orbit density correction is given as: 




Due to the spin degeneracy, various correction densities 
should be multiplied by a factor of 2. Note that there is no 
zeroth-order spin-orbit density. Since the liner term of (p x a) 
doesn’t contribute in the integration, the final expression of 
J 2 (r) can be obtained by considering only the Gi_so’s contri¬ 
bution. 
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( 21 ) 


Up to now, we have not explicitly considered the density- 
dependent effective mass. The density-dependent effective 
mass is involved in the central part of the single-particle 
Hamiltonian, 

hc = ^f{r)p'^ + U{r) with m*(r) = l//(r). (16) 

In fact, the expressions of the 2'''^-order Thomas-Fermi do not 
need to be modified. However, three additional expansion 
terms appear compared to the cases without effective mass: 

h1i^h = fpM^f + 2fVM-^p^iVM (17a) 

h{Mh)^ = -l/p4(v/)2 + /(vc/)2 + i/Vv2/ 

+ i/Vv2t7 + i/p2(V/)(V[/) 

(17b) 

hC^h){tA) = i/p2(v/)(VA) + /(VC/)(VA) (17c) 

Finally, the normal density (and pairing density) in real 
space can be obtained by integrals over quasiparticle energies 
E or single-particle energies e up to a certain cutoff 


III. THE LIMIT OF ZERO PAIRING GAP 


In the limit of the zero pairing gap, the derived general¬ 
ized 2'''^-order Thomas-Fermi approximation for the normal 
state should be recovered. For the normal density and pairing 
density, this has been demonstrated in Ref. In this work, 
we have to examine the additional terms of spin-orbit con¬ 
tributions and the spin-orbit density J(r) in the limit of zero 
pairing gap. For the spin-orbit contribution to normal density: 


lim 
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AE^ 



( 22 ) 


The resulting density expression is equivalent to the 2"*^-order 
Thomas-Fermi approximation of the non-superfluid state ll^ . 
Obviously, the resulting pairing density due to the spin-orbit 
contribution becomes zero in the limit of the zero pairing gap. 
Similarly, the spin-orbit density in the limit can be obtained. 
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which also agrees with the expression of the normal state 
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IV. IMPLEMENTATION AND CALCULATIONS 

To demonstrate applications of the 2“'^-order Thomas- 
(18) Fermi corrections, we perform self-consistent calculations of 

238 

a finite deformed nucleus U. For the particle-hole inter¬ 
action channel, the often used Skyrme force SLy4 El is 
adopted. For the pairing channel, the volume pairing, i.e., a 
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FIG. 1. (Color online) The 0**'-order (dotted line) and 2“‘^-order 
Thomas-Fermi contributions (0*^-order plus 2"'^-order ) to the neu¬ 
tron pairing density from 55 to 65 MeV, as well as the correspond¬ 
ing coordinate-space HFB solutions (as labeled by ‘Box’). The 2“'^- 
order Thomas-Fermi contribution without explicitly considering the 
effective mass is also shown, as labeled by ‘2TF-(w/o f)’. 



R(fm) 

FIG. 2. (Color online) The 2“‘^-order Thomas-Fermi terms corre¬ 
sponding to the pairing density in Fig[T] The terms denoted by 9i as 
listed in Eq. 01. The spin-orbit contribution (dashed line) and the 
summed correction (dotted line) are shown. 


delta interaction for the pairing force is adopted with a reason¬ 
able pairing strength of 200 MeV fm“^. The Skyrme density 
functional contains the density-dependent effective mass and 
the spin-orbit potential, providing an ideal testing ground for 
the generalized 2''‘^-order Thomas-Fermi approximation. 


A. Non-self-consistent Calculations 

The nuclear HFB solutions have deep-hole states S 
which correspond to deep-bound states in Hartree-Fockn-BCS 
solutions. In the HFB approach, deep-hole states are narrow 
quasiparticle resonances with large quasiparticle energies 
and large occupation numbers v\. These states contain impor¬ 



R (fm) 

FIG. 3. (Color online) Similar to Fig[T] but for the 0*''-order and 
2"'^-order Thomas-Fermi contributions to the neutron normal den¬ 
sity. The normal density from self-consistent coordinate-space HFB 
solutions is labeled by ‘Box’. 



FIG. 4. (Color online) Similar to Fig|^ but for the 2"‘^-order contri¬ 
bution terms to the normal density as listed in Ea.lll3b. 


tant shell effects and can not be well described by the Thomas- 
Fermi approximation, since the latter does not account for 
shell effects. Therefore, we first aim to compare the con¬ 
tributions of various partial densities corresponding to high 
quasiparticle energies. Our final goal is to combine the HFB 
solutions at low-energy quasiparticle states and 2'''^-order cor¬ 
rections in the quasiparticle high-energy region. In such a hy¬ 
brid way, the computational costs for large superfluid Fermi 
systems can be remarkably reduced. Indeed, the high-energy 
states behave like quasi-continuum states and are distinctly 
different from low-energy states. One can suppose that the 
former can be well approximated by the Thomas-Fermi ex¬ 
pressions. 

First, self-consistent HFB calculations are performed for 
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U, then the zeroth-order and second-order Thomas-Fermi 
approximations are investigated by performing one iteration 
with densities from full HFB solutions. The HFB calcula¬ 
tions are performed in the cylindrical coordinate space with 
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FIG. 5. (Color online) Similar to Fig[T] but for the 0*''-order and 2“'^- 
order Thomas-Fermi contributions to the neutron kinetic density. 



FIG. 6. (Color online) Similar to Fig[T] but for the 2“‘^-order 
Thomas-Fermi contribution (0*^-order doesn’t contribute) to the di¬ 
vergence of the neutron spin-orbit density V ■ J. 


the HFB-AX solver that uses B-spline techniques. For 
calculations employing large box sizes and small lattice spac- 
ings, the discretized continuum spectra would be very dense, 
providing good resolutions. The pairing densities p{r, z) are 
plotted along the diagonal coordinate R — Vr^~+~^\{r=z) 
rather than the axes, to avoid numerical errors at boundaries. 
We compare various pairing densities in Fig[T] In this work, 
our discussions are restricted to neutrons since pairing is ah- 

238 

sent in protons in U. The pairing densities of the full self- 
consistent HFB calculation are summed for quasiparticle ener¬ 
gies ranging from 55 to 65 MeV. This is the high energy win¬ 
dow for which we want to study the Thomas-Fermi approx¬ 
imation, in discarding the influence of deep-hole states. The 
upper cutoff is taken as 65 MeV in all our results. The con¬ 
verged FIFB densities are then inserted into the Thomas-Fermi 
expressions for one iteration and the zeroth and second order 
Thomas-Fermi densities are also summed (integrated) over 
the same energy interval. In this way the quantities shown 



FIG. 7. (Color online) (a) The 2“‘^-order Thomas-Fermi terms to 
the neutron normal density with an energy cutoff from 25 MeV to 65 
MeV; (b) the 2"‘^-order terms to the pairing density; (c) the 2“‘^-order 
terms to the pairing density with a pairing strength 10 times larger. 


in Fig[T]are perturbative. It can be seen that the second-order 
contribution is very small but can slightly improve the pairing 
density distribution. Indeed, as shown in Eg.lfTS]) and Ea. (fT4l) . 
the second-order corrections involve terms of order 1/E^ and 
1 / E"^, which are suppressed in the high-energy region. With 
the derivative terms of effective mass, the second-order cor¬ 
rections are only slightly modified in the high energy region. 

In Figl2l the second-order contributions from different 
terms in Eq. (fT4l i are displayed corresponding to EiglT] As we 
can see, the spin-orbit effect has non-negligible contributions. 
The term of h A is dominating in the second-order pair- 
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FIG. 8. (Color online) The neutron pairing density corresponding to 
FigOc) plotted on a logarithmic scale to show surface distributions, 
compared to the coordinate-space HFB solutions. 


FIG. 9. (Color online) The hybrid HFB calculations of neutron pair¬ 
ing density in o*^-order and 2“‘^-order Thomas-Fermi ap¬ 

proximations within an quasiparticle energy interval from 25 MeV to 
65 MeV, compared to the HFB-AX calculations. 


ing density. In the coarse graining treatment of Bogoliubov- 
de Gennes equations, only this term has been included in the 
pairing equations iflOl] . Based on our full analysis, we see 
that such a coarse graining treatment is reasonable, in particu¬ 
lar at the surface. The terms of h{ A h)^ and h A'^h are non- 
negligible but the two almost cancel each other in the high 
energy region. 

Similar to FiglT] the contributions to normal densities are 
displayed in Figl3 Compared to FiglT] the high-energy con¬ 
tributions to the normal density are much smaller (by 2 orders 
of magnitude) than to the pairing density. The second-order 
correction can improve the density distribution at the surface. 
However, the second-order correction leads to enhanced oscil¬ 
lations in the inner part of the density distribution, due to non- 
self-consistent calculations. In FigHj we see that the terms 

^ _ y _ y ^ 

of A( A /i)^ and h AA are dominating but the two almost 

^ ^ _ y ^ _ y 

cancel each other. The terms of h A^h and h{ A Kf' are also 
canceling each other. The term of h{ A A)^ is negligible. In 
contrast to the pairing density, it is hard to say which term 
plays a major role in the normal density. 

In Fig|2 the corrections to kinetic densities are displayed. 
Again, the second-order correction improves the description 
at surfaces compared to the zero-order correction. In addition, 
the second-order correction also improves the kinetic density 
in the inner part. In Fig|6l the second-order correction to the 
divergence of the spin-orbit density, V • J, is shown. Note that 
there is no zero-order correction to the spin-orbit density. The 
obtained V • J roughly agrees with the HFB solutions. It is 
important to obtain the spin-orbit density to do self-consistent 
semiclassical calculations. 

To study the cutoff dependence, the second-order Thomas- 
Fermi correction terms to the densities with an enlarged en¬ 
ergy interval from 25 MeV to 65 MeV are shown in Fig|7| 
In FiglTja), the corrections to normal densities are displayed. 
Compared to FigH) the resulting densities are increased by 
a fqcjt or of 20. It is obvious that the terms of h A^h and 
h{ A h)^ do not cancel each other anymore. The same holds 


^ _ y ^ _ y ^ _ y 

for the terms of A( A /i)^ and h A ^A. The term h{ A A)^ is 
still negligible. In Fig|2]b), the corrections to pairing densi¬ 
ties are shown. Compared to Fig|2l the resulting pairing den¬ 
sities are increased by a factor of 10. In this case, we see 

^ _ y 

that the total correction do not follow the term ft, A ^A any¬ 
more which was dominating in the high energy window. The 
terms of ft( A ft)^ and ft A ^ft do not cancel each other. Or, in 
other words, the coarse graining treatment is not applicable in 
the low-energy region. The terms of second-order have to be 
fully taken into account for low-energy quasiparticle states. 

In fact, the pairing in nuclei is rather weak, compared to 
cold atomic systems, for instance, close to the unitary limit. 
Therefore it is instructive to study for what happens if we in¬ 
crease artificially the pairing. In FiglT^c), the resulting cor¬ 
rections to the pairing densities are obtained with a 10-times 
larger pairi^ strength compared to FiglT^b). In this case, the 
term of A( A ft)^ acquires a role which is absent in Fig|2tb). 

In FigO the pairing density contribution corresponding to 
FiglTJc) are shown on a log scale. It can be seen that the 
second-order correction agrees exactly with the asympotics 
of coordinate-space HFB solutions at large distances. While 
the zeroth-order correction underestimates the pairing density 
at the nuclear surface by 10%. These examples demonstrate 
clearly the advantages of the full second-order Thomas-Fermi 
approximation for quasiparticle energy intervals with a lower 
edge reaching into the low energy domain (e.g., 25 MeV). 


B. Hybrid HFB Calculations 

One of our main motivations to implement the hybrid HFB 
calculations with the second-order Thomas-Fermi approxima¬ 
tion is to reduce the computing costs for large systems. In 
our previous work, we have tested the hybrid HFB strategy 
with the zeroth-order Thomas-Fermi approximation ra . In 
the hybrid strategy, the high energy deep-hole states and con- 
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FIG. 10. (Color online) The hybrid HFB calculations of binding en¬ 
ergies of with 0*’'-order and 2"‘^-order Thomas-Fermi approx¬ 
imations at different cutoffs, compared to the HFB-AX calculations. 

tinuum states are separately treated. In the present work, the 
high-energy continuum is treated in the second-order Thomas- 
Fermi approximation. The deep-hole states, which are narrow 
quasiparticle resonances, can be described separately. Firstly, 
we diagonalize the single-particle Hamiltonian to obtain the 
deep-bound Hartree-Fock wavefunctions and single-particle 
energies. 

{h-X)vf^{r)=e,vf^{r) (24a) 


{h- \ + ei)uf^{r) = A(r)t;f ^(r) (24b) 

Next, one natural way is to renormalize the deep-hole wave- 
functions with the BCS approximation, as in Ref. fH. We 
solve the Ea. (l24bb to describe the asympotics of scattering 
components in stead of BCS. Another more elaborate way 
is to solve the HFB equation perturbatively, as described in 
Ref. m. 

In the hybrid strategy, we introduce an energy interval (or 
window) of quasiparticle energies from a lower cutoff to 65 
MeV. Below the lower cutoff, the box discretization is used 
to solve the HFB equation using the HFB-AX solver. Within 
the window, the zeroth-order or second-order Thomas-Fermi 
approximation for the continuum and the BCS approximation 
for deep-hole states are used. The quantal effects are taken 
into account by the standard HFB solutions below lower 
cutoffs and the deep-hole states. The hybrid HFB calculations 
are implemented self-consistently to keep the conservation 
of particle numbers. With iterations in terms of various 
densities, i.e., the normal density p{r), the kinetic density 
rir), the spin-orbit density J(r) (does not appear in the 
zeroth-order), the pairing density p{r), the self-consistent 
hybrid Skyrme HFB calculations can be realized. The quantal 
and semiclassical are coupled via the summed densities. In 
our test calculations of we use the SLy4 force and 

the volume pairing interaction, as we used in the non-self- 
consistent calculations. In Fig|9] the total pairing density 


prohles from hybrid calculations are shown. It can be seen 
that the agreements between hybrid and full HFB calculations 
are quite good. 

Concerning the treatment of the deep hole states, we should 
note that the BCS underestimates the pairing correlation com¬ 
pared to the HFB approach, due to the absent of continuum 
coupling. For example, to reproduce the neutron pairing gap 
of 1.245 MeV in ^^°Sn with SLy4 and the volume pairing, 
the pairing strengthes in BCS and HFB have to be adjusted 
to 283 and 187 MeV fm“^, respectively. Hence we slightly 
increase the pairing in the BCS treatment of deep-hole states 
with one global factor. As a result the average pairing gaps at 
all cutoffs are close to the full HFB result. For systems such 
as trapped Fermi gases, there are no deep-hole states and the 
hybrid calculations will be simpler. 

The results of hybrid HFB with different lower cutoffs of 
the energy windows are displayed in Fig(T0] In Fig(T0l the 
deviations in total binding energies between the hybrid HFB 
and the coordinate-space HFB solutions increase as the lower 
cutoff decreases. Due to the self-consistency, the second- 
order approximation is obviously better than the zero-order 
Thomas-Fermi, in particular at low-energy cutoffs. As we can 
see, with the cutoff energy at 25 MeV, the deviation is ~0.5 
MeV over the total binding energy of 1790.52 MeV, which is 
quite satisfactory. Note that the deviations are not only from 
the approximate treatment of continuum but also from the 
treatment of deep-hole states. In the low-energy region, quasi¬ 
particle resonances can acquire considerable widths, leading 
to ambiguous contaminations. In addition, the numerical ac¬ 
curacy of derivatives is important in the second-order calcula¬ 
tions. It will be improved by the multi-wavelet techniques ifisll 
in the future. The full 3D coordinate-space HFB calculations 
are very expensive even with the efficient multi-wavelet tech¬ 
niques id. With the lower cutoff of 25 MeV, the number of 
eigen-functions to be solved would be reduced by half. On the 
other hand the extra numerical cost to include the second order 
Thomas-Fermi term with respect to take only the zeroth-order, 
is almost negligible. In this case, the resulting computational 
cost can be reduced at least by one order of magnitude and the 
desired accuracy is still retained. 


V. SUMMARY 

In summary, we extended the second-order Thomas-Fermi 
approximation of the Hartree-Fock-Bogoliubov solutions for 
superfluid systems by including the effective mass and the 
spin-orbit potential, which are essential for full calculations. 
In particular, the spin-orbit contribution is important because 
it is absent in the zeroth-order Thomas-Fermi approximation. 
The expressions of the new terms have been checked in the 
zero-pairing limit. The full second-order terms have been ex¬ 
amined numerically in a perturbative manner in comparison 
with self-consistent coordinate-space HFB solutions. In gen¬ 
eral, the second-order corrections can improve various density 
distributions at surfaces compared to the zeroth-order correc¬ 
tions. The signihcance of including full superfluid second- 
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order corrections definitely increases as the lower edge of the 
energy interval in which the Thomas-Fermi approximation 
acts goes down. Among the second-order correction terms, 
the pairing density contribution is relatively the most impor¬ 
tant one. Furthermore, we performed fully self-consistent hy¬ 
brid HFB calculations with the second-order Thomas-Fermi 
for the continuum and satisfactory results have been obtained. 
This will be particulary useful for the 3D coordinate-space 
HFB calculations which are computationally very expensive. 
The second-order superfluid Thomas-Fermi method can be 
further extended by including the temperature dependence and 
rotational or vector fields. More interesting applications of 
the generalized second-order Thomas-Fermi approximation 
are intended for large complex inhomogeneous superfluid sys¬ 
tems such as cold atomic condensates and neutron star crusts, 
in which the pairing fields are very large compared to hnite 
nuclei. 


Eq. (fTSl l are listed here: 


Vi = 

m ■- 

V3 -- 

m = 
Vs = 
V6 = 
V7 ■- 
V8 = 
V9 = 
Vio 


3hA2 

~32E^ 

{2h^-A^)A 

16E^ 

2/iA2 - 
32^5 

4^2 A2 - A4 

16^7 

(2^2 - 3A2)/iA 

8^7 

5h2A2 

16^ 

2hA^ - 3h^A 

16^7 

3^2 A2 -h^-A^ 

8^7 

2hA^ - 3hA^ 
16F;7 
+ hA2 
^ 16F;7 


(A.l) 


The expansion coefficients of the second-order super¬ 
fluid Thomas-Fermi approximation to the pairing density in 
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Appendix: Second-order Expansion Terms 


0i 

02 

03 


05 

06 


07 


08 


09 


(A2 - 2h^)A 
32^5 

/l(h2 - 2A2) 
16E^ 

3/i2A 
~ 32E^ 

(3A2 - 2h^)hA 

16^7 

hi + A^- 3^2 A2 

8^7 

(2A2 - 3h^)hA 

16^7 

5/i2A2 

16£;7 

(2A2 - 3h2)hA 

32^7 

h* - 4^2 A2 

16^7 


^10 = 


A 

16E^ 


(A.2) 


The non-zero derivative terms employed in Eq. ( fT3] ) and 
Eq. (HI: 


hA^A = fV^A (A.3a) 


The expansion coefficients of the second-order super¬ 
fluid Thomas-Eermi approximation to the normal density in 


+ 2fV^U - |p2(V/)2 (A.3b) 
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{A^hf = ^fp\VAf (A.3c) 


Ai^hf = ^fp\V^A) (A.3f) 


h{^h)^ = -l/p4(v/)2 + /(vc/)2 + l/VvV 
+ i/Vv2[/ + l/p2(v/)(VC/) 

(A.3d) 

hi^h)C^A) = i/p2(v/)(VA) + /(VC/)(VA) (A.3e) 
o 


h{ A A)2 = /(VA)2 (A.3g) 
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